Gap solitons in quasiperiodic optical lattices 
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Families of solitons in one- and two-dimensional (ID and 2D) Gross-Pitaevskii equations with the 
repulsive nonlinearity and a potential of the quasicrystallic type are constructed (in the 2D case, 
the potential corresponds to a five- fold optical lattice). Stable ID solitons in the weak potential are 
explicitly found in three bandgaps. These solitons are mobile, and they collide elastically. Many 
species of tightly bound ID solitons are found in the strong potential, both stable and unstable 
(unstable ones transform themselves into asymmetric breathers). In the 2D model, families of both 
fundamental and vortical solitons are found and are shown to be stable. 

PACS numbers: 42.65.Tg, 03.75.Lm, 61.44.Br, 42.70.Qs 

I. INTRODUCTION 

Solitons in Bose-Einstein condensates (BECs) are currently a subject of intensive theoretical and experimental 
studies. Solitons supported by weak attractive interactions between atoms were created in the condensate of 7 Li 
trapped in strongly elongated (nearly one-dimensional, ID) traps 1] (although the actual shape of the solitons 
was actually nearly three-dimensional, rather than nearly-lD; the latter feature was observed in a salient form in the 
recent experiment 2], where solitons were created as a residual pattern after dynamical collapse took place in the 85 Rb 
condensate). The stability of the solitons, with a given number of atoms, against collapse is secured, in these settings, 
by a combination of the trap's geometry and small absolute value of the negative scattering length characterizing the 
attraction between atoms (~ 0.1 nm, in the case of 7 Li, but up to ~ 2 nm, in the condensate of 85 Rb). A very accurate 
description of the solitons is provided by the mean-field approximation based on the Gross-Pitaevskii equation (GPE) 



Positive scattering length, which corresponds to repulsive interactions, is more generic in BEC (in the above- 
mentioned experiments, the negative scattering length was actually artificially induced by means of the Feshbach 
resonance, both in 7 Li and 85 Rb). In a repulsive condensate, gap solitons (GSs) may be created as a result of 
the interplay of the self-defocusing nonlinearity and periodic potential induced by an optical lattice (OL, i.e., an 
interference pattern created by counterpropagating laser beams illuminating the condensate) 0, || . GSs emerge in 
bandgaps of the system's linear spectrum, since the combination of a negative effective mass, appearing in a part of 
bands adjacent to the gaps, with the repulsive interaction is exactly what is needed to create a soliton. Theoretical 
models for GSs in BEC were reviewed in Ref. |g, and rigorous stability analysis for them in the ID case was developed 
in Refs. [jj. Creation of a GS in the experiment was reported in Ref. |8|, in 87 Rb condensate placed in a quasi-lD 
trap supplemented by a longitudinal OL; the soliton contained a few hundred atoms. In a subsequent experiment, 
large-size confined states with much larger numbers of atoms were discovered in a stronger OL 9] ; an explanation to 
this observation was recently proposed 10], which, essentially, treats the extended state as a segment of a nonlinear 
Bloch wave, bounded by two fronts (domain walls) which are sustained by the strong OL (a similar wall between filled 
and empty domains was predicted in BEC with self-attraction in Ref. pj|). 

Stable GSs were also predicted in 2D and 3D settings 0, including 2D solitons with embedded vorticity [2, EH- In 
nonlinear optics, similar 2D spatial solitons 14] and 2D localized vortices 15] were predicted in photonic crystals and 
photonic-crystal fibers, as well as in periodic photonic lattices induced by perpendicular laser beams in photorefractive 
materials Il6|. In media of the latter type, both fundamental [l7| and vortical 2D solitons have been created in the 
experiment [I3 • A difference from the self-repulsive BEC is the self- focusing character of the nonlinearity in nonlinear 
optics (which is cubic in photonic crystals, and saturating in photorefractive media); for this reason, the optical solitons 
usually belong to the semi-infinite bandgap, where the effective mass is always positive (although vortex solitons in a 
finite gap were also created in a photorefractive medium [19^. 

The GPE with self-attraction and periodic (optical-lattice) potential gives rise to 2D solitons and vortices [2(J 
similar to those found in the above-mentioned optical models. Moreover, it has been recently demonstrated that both 
the GPE with the self-focusing cubic term and OL potential, and its counterpart with the saturable nonlinearity, that 
appertains to photorefractives, give rise to stable higher-order localized vortices (which are built as rings of unitary 
vortices) and supervortices (similar rings with global vorticity imprinted onto them) [2]] . 
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Quasiperiodic OLs can be easily created too: in the ID case, as a superposition of two sublattices with incommen- 
surate periods, and in the 2D case, as a combination of N = 5 or N > 7 quasi- ID sublattices with wave vectors 
of equal lengths, which make equal angles 2tt/N between themselves. In particular, the 2D lattice with N = 5 is 
known as the Penrose tiling (PT). The bandgap spectrum of 2D photonic crystals of the PT type has been studied 
in detail with a conclusion that true (omnidirectional) bandgaps may be supported by the PT. 

The use of quasiperiodic lattices, in one and two dimensions alike, offers new degrees of freedom that allow one 
to engineer desirable bandgaps in the spectrum and, this way, design soliton families in nonlinear media. Several 
theoretical works sought for solitons in models combining quasiperiodic lattice potentials and cubic self -focusing. In 
an early work 23], solitons were not found in a ID quasiperiodic model; however, they were later discovered in a 
"deterministic aperiodic" discrete nonlinear Schrodinger (NLS) equation, which may be regarded as a limit case of 
the ID model with a very strong quasiperiodic potential [24|. Localized and delocalized solutions were also studied 
in the ID model with a superlattice (i.e., an ordinary OL subjected to an additional modulation, which is different 
from a quasicrystal) |2£|. Recently, 2D solitons (with zero vorticity) were numerically constructed in a model of a 
photonic crystal made of a self-focusing material, with N = 12, in terms of the above definition 26]. 

The objective of this work is to find GSs supported by the interplay of the cubic repulsive nonlinearity and quasiperi- 
odic lattice potentials. We will report systematic results for fundamental solitons in the ID and 2D models (the latter 
one will be elaborated for the PT lattice). In the ID case, loosely bound GSs in the weak lattice are mobile. Tightly 
bound GSs in the strong lattice are found in many modifications, some stable and some not. Stable vortex solitons 
in the PT potential will be demonstrated too. 

II. ONE-DIMENSIONAL SOLITONS 

In the mean-field approximation, the evolution of the single-atom function <j) obeys the GPE with the repulsive 
nonlinear term and quasiperiodic potential U (x) . In the normalized form, the equation is 

As said above, the ID potential is a combination of two incommensurate spatial harmonics with equal amplitudes e, 

U(x) = -e[cos (tt(x - L/2)) + cos (qir(x - L/2))), (2) 

where x = L/2 is the central point of a large(L ^> 1) trapping domain, < x < L, the period of the first sublattice is 
normalized to be 2, and q is an irrational number. Below, we display results for q = (>/5 + l)/2 w 1.62. 

As is known, an exact bandgap spectrum of the linear Schrodinger equation with a quasi-periodic potential is 
fractal. Without the aim to display the spectrum in full detail, in Fig. 1 we show its part which is relevant to the 
quest for gap solitons. Bands of values of chemical potential \i corresponding to families of quasi-Bloch states, 

4>(x,t) = e-^xix), (3) 

are covered by vertical segments at fixed values of OL strength e. In panel 1(a), one can discern four finite gaps 
(discontinuities between the bands, the lowest gap being extremely narrow) and the underlying semi-infinite gap (one 
extending to \i — > — oo). A result of the standard perturbation theory is that these gaps start (at e = 0) at fi = tt 2 /2 
and \i = (irq) 2 /2. Panel 1(b) details the results in a narrower range of \i but for an extended interval of e. One can 
observe, in particular, that the lowest band in Fig. 1(a) splits into two bands at e = 0.8, and three bands at e = 1 in 
Fig. 1(b). With the increase of e, additional bandgaps open up on a finer scale (accurate computation of the spectrum 
becomes rather difficult for e > 2). 

As the effective mass is negative near the top of (quasi-)Bloch bands, GSs are expected to exist in adjacent 
lower parts of the gaps. Assuming real \{ x ) an d varying \i in Eq. (|3J), we numerically searched for a family of 
solutions of the corresponding stationary version of Eq. JI}, x"/2 = x 3 + U ( x )x ~ A% satisfying boundary conditions 
x(x = L/2) = A, x'( x = L/2) = x( x = L) = 0. Figure 2 displays three typical examples of GSs with equal amplitudes, 
x(x = L/2) = A = 0.3, that were found in three finite bandgaps for e = 1 at values of \i marked by crosses in Fig. 
1(b) [as usual (in the GPE with self-repulsion), no solitons in the are found in the semi-infinite gap]. A blow-up of 
the GS in Fig. 2(b) is shown in Fig. 2(d), to demonstrate that, while the envelope of the solution is that of a solitary 
wave, the carrier wave function is nearly quasiperiodic. 

Families of the GS solutions are characterized by the dependence of /i on the soliton's norm, N = x 2 ( x )dx. 
For GS families in the three lowest finite bandgaps, whose examples were displayed in Fig. 2, the N(ja) dependences 
are plotted in Fig. 3. Naturally, N — > corresponds to \i approaching the border of the quasi-Bloch band located 
under the respective gap. 
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FIG. 1: A part of the bandgap structure in the linear version of Eq. (IJ. Crosses in (b) indicate values of \i for three GSs 
shown in Fig. 2. 
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FIG. 2: Examples of stable gap solitons in the ID model for e = 1, with the chemical potential and norm fi — 0.658, N — 0.104 
(a), fj, = 0.237, N = 0.763 (b), and ji = 0.0196, N = 0.935 (c). Panel (d) is a blow-up of (b) in the region of 150 < x < 190. 



The numerically found solution branches do not cover the range of \i corresponding to the entire bandgap in each 
case. However, the effective termination of the branches inside the gaps seems to be a numerical problem, rather than 
a true feature of the model: it is difficult to secure convergence of the numerical solutions for very large values of N. 

Stability analysis of the GSs was performed by means of direct simulations, using the split-step method with 4096 
Fourier modes. This way, it was concluded that all the three families displayed in Fig. 3 are completely stable. 

GSs obtained above can be readily set in motion the same way as it was done in Ref. [27|, i.e., multiplying a 
stationary solution by exp(zfc#), with k not too large. Examples are displayed in Fig. 4, where the "shove factor" 
and the resulting average velocity of the moving soliton are k = 0.03, v = —0.104 (a); k = 0.015, v = —0.117 (b); and 
k = 0.005, v = -0.03 (c). Accordingly, the effective GS mass M = k/v is -0.288 (a), -0.128 (b), and -0.17 (c). The 
mass is negative, as it should be for GSs 27]. On the other hand, moving solitons feature conspicuous radiation losses 
if they are set in motion by shove k exceeding a certain critical value, which is k cr « 0.05, k « 0.017, and k ~ 0.008 
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FIG. 3: The chemical potential vs. norm for gap-soliton families in three lowest finite bandgaps (a,b,c) of the ID model with 
e=l. The dashed curves are added as guides to the eye. 
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FIG. 4: Examples of stable moving gap solitons in three finite bandgaps (a,b,c) of the ID model. 
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FIG. 5: Head-on collision of two gap solitons, which are obtained from ones displayed in Fig. 2(a) by shoving them with 
exp(±ikx), k = 0.02. 



for the solitons from Figs. 2(a), 2(b), and 2(c) respectively, although the transition to the lossy motion regime is not 
very sharp. 

Following the approach based on the separation of rapidly and slowly varying functions, one can approximate GS 
solutions by 



(4) 



where fn(x) is a linear quasi-Bloch function pertaining to chemical potential ji, and &(x,t) is a slowly varying 
amplitude, for which an averaged NLS equation can be derived as in Ref. |27j. 



<9$ 



1 <9 2 $ 
' 2M dx 2 



G|$| 2 $. 



(5) 



Here M < is the above-mentioned effective mass, and G > is an effective nonlinearity coefficient, which can be 
found by matching obvious soliton solutions to Eq. JJJ, $ = Asech [Ay/— MGx) (A is an arbitrary amplitude), to the 
envelope of the numerically found GSs. As a result, for the three quiescent GSs displayed in Fig. 2 we find G = 0.68 
(a), G — 0.41 (b), and G = 0.13 (c). The availability of stable moving solitons suggests to consider collisions between 
them. As predicted by the averaged equation J5J), moving GSs emerge unscathed from collisions, see an example in 
Fig. 5. 

The GSs displayed above are loosely bound solutions, in terms of Ref. [27|, as they are generated by the GPE with 
a relatively small strength (e) of the pinning potential. For larger e, bands between the gaps in the model's spectrum 
become nearly invisible, and the character of GS solutions drastically alters, with a large variety of tightly bound 
solitons found in the broad bandgaps. Figure 6(a) displays four different kinds of the solutions (for e = 5) with equal 
amplitudes, \4>(x = L/2)\ = 2. In the first potential well around the central point, all the profiles are almost identical, 
but the next pair of peaks appear at markedly different positions, and may even have different signs. The form of the 
quasiperiodic potential U(x) is displayed below the GSs. Naturally, local maxima of \(j)\ are found near minima of the 
potential. In fact, more GS species can be found, in addition to the four ones shown here. Figure 6(b) quantifies four 
GS families, whose representatives are displayed in Fig. 6(a), in terms of the fi(N) dependence. The curves overlap 
at sufficiently small N, as in this limit each GS species virtually reduces to the single peak in the first potential well; 
however, the species become very different at larger A/", corresponding to the different GS shapes displayed in Fig. 
6(a). For example, branch 4 starts to deviate from the other three ones near N = 0.5, when negative peaks begin to 
develop near x = 16 and 24. 
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FIG. 6: (a) Four types of stationary gap solitons found in the ID model with a strong optical lattice, e = 5. The chemical 
potential \i and the norm N are, respectively, /x = —1.261, N = 1.853 (1), \i = —1.088, N = 1.266 (2), fi = -1.214, N = 1.697 
(3), \i — —1.2415, N = 2.854 (4). The form of the potential is shown in the bottom. Panel (b) displays curves fi(N) for the 
four species of the gap solitons. Solution branches shown in (b) terminate at points where the solitons quickly develop the tail 
structures. 
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FIG. 7: (a) Evolution of \<f>\ for an unstable gap soliton labeled by 3 in Fig. 6(a). (b) Evolution of \<f>(L/2)\ in the gap soliton. 



Direct simulations demonstrate that the GS species labeled as 1, 2 and 4 in Fig. 6 are stable, while the family labeled 
by 3 is unstable. Further, Fig. 7(a) displays a typical example of the evolution of an unstable GS. The instability breaks 
the soliton's symmetry and makes it a breather. Although the breathers's amplitude features complex oscillations, 
see Fig. 7(b), the breather does not decay, maintaining its localized shape. 



III. TWO-DIMENSIONAL SOLITONS 



In the normalized form, the 2D version of the GPE with the PT (Penrose tiling) potential is 

5 



dt 2 



101' 



cos (k< n > • r) 



(6) 



with |fci n \fc^| = 7r {cos (27r(n — l)/5) , sin (27r(n — l)/5)}. GS solutions in two dimensions were constructed by 

means of a combination of the known method of the integration in imaginary time [28;) and subsequent real-time 
simulations. To this end, a substitution was first made, <fi(x,y,t) = e~ ZfJjt $>(x, y, — ir), which transforms Eq. © into 
an nonlinear diffusion equation, 



d<$> 1 O 9n 

— = -V 2 $ + (/i -$ 2 )$ + £ 



cos (k (n) • r) 



(7) 



This equation was solved numerically by dint of the split-step Fourier method with 256 x 256 modes, in a domain of 
the size of 50 x 50. It was observed that the norm of the solution originally decreased, and then began to increase. 
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FIG. 8: Contour plots of \<p(x,y)\ for three examples of stable gap solitons in the 2D model with a strong Penrose-tiling optical 
lattice, £ = 5. The solitons are generated by initial conditions 0, with, respectively, \i — —1.56 (a), \i — —3.33 (b), and 
fi = -0.41 (c). The norms of the three solitons are N a = 147, N b = 103, and iV c = 268. 




FIG. 9: Time dependence of the peak amplitudes of three localized solutions shown in Fig. 8. 



The imaginary-time integration was switched back into simulation of the GPE in real time when the norm attained 
its minimum. 

First, we present typical examples of stable fundamental (zero-vorticity) 2D solitons, for \i — —1.56, —3.33, and 
—0.41. They were generated by initial conditions [in Eq. Q] with, respectively, 

f exp(-0.5(x 2 + ^ 2 )) ] 5 f v 
$o(*,2/)=< 1.5exp(-O.150r 2 + 2/ 2 )) } V cos -k^ • r , (8) 
[ 2exp(-0.05(x 2 + y 2 )) J n=i ^ ' 

the last multiplier being a half-harmonic counterpart of the PT potential in Eq. (|BJ), which is expected to parametrically 
couple to the potential. 

Figure 8 displays contour plots of \<fi(x,y)\ in final states produced by the numerical integration. To illustrate 
the proximity of the solutions to truly stationary ones, and their stability, in Fig. 9 we show the asymptotic time 
dependence of the peak amplitudes of the three solitons at (x,y) = (L/2,L/2) = (25,25). These gap solitons are 
found in a strong optical lattice with e = 5, for which bands separating the spectral gaps are extremely narrow. 

Systematic results for fundamental 2D solitons were generated by using initial conditions (|BJ) and varying the 
respective chemical potential \i. Figure 10 shows /jl(N) dependences for three families of thus generated GSs, where 
N = J J \(j)(x,y)\ 2 dxdy is the usual 2D norm. At relatively small values of TV, the dependences completely overlap, 
and they slightly split up at very large iV, corresponding to \i — > —0. They are somewhat (but not quite) similar to 
plots for tightly bound GSs in the ID model, cf. Fig. 6(b). There are no loosely bound stable GSs in two dimensions, 
because the averaged two-dimensional NLS equation [a 2D counterpart of Eq. ®] does not have stable soliton 
solutions. 

Lastly, stable 2D solitons with embedded vorticity S can be found too. They were generated by adding a factor, 
(x 2 + y 2 ^ s ^ 2 ex.-p(iS0), which distinguishes vortex states in uniform media, to initial ansdtze (|BJ). A typical example 
of such a stable vortical gap soliton in presented Fig. 11. While the field pattern of |(/>(x, y)\, displayed in panel 11(a), 
is stationary, panel 11(b) is actually a snapshot, as the pattern of Re <f>(x,y) shown in this panel rotates clockwise, 
at the angular velocity of uj = n/ S = — 3. 
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FIG. 10: Chemical potential fi of two-dimensional fundamental solitons vs. N, in the model with e = 5. Rhombuses, crosses, 
and squares represent soliton families generated, respectively, by three initial configurations in Eq. (0 (typical examples of 
solitons belonging to the three families are given in Fig. 8). 
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FIG. 11: A typical example of a stable vortical gap soliton with S = —1, found in the model with s = 5. (a) and (b): 
Contour plots of \cf>(x,y)\ and Re <j)(x,y), the latter shown only in the region with Re <p(x,y) > 0. (c): A sequence of central 
cross- sect ions, \<j)(x — L/2,y)\, taken at different moments of time, which demonstrate the vanishing of the field at the center, 
as must be in the vortex, and its stability. The norm and chemical potential of the solution are N = 254 and fi = — 3. 



IV. CONCLUSION 



We have constructed families of ID and 2D gap solitons in the model combining the self-defocusing nonlinearity and 
quasiperiodic lattice potentials. Loosely bound ID solitons supported by a weak lattice were found in three gaps; they 
are mobile, and collide elastically. The strong ID lattice supports a variety of species of tightly bound solitons. Most 
of them are stable, unstable ones giving rise to robust asymmetric breathers. 2D gap solitons, both fundamental and 
vortical ones, are stable too, and they may only have a tightly-bound shape. The predicted solitons can be created by 
means of available experimental techniques in a Bose-Einstein condensate with a positive scattering length, loading 
it into an optical lattice induced by a superposition of several laser beams. 
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